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Abstract 

Light propagation in materials with microscopic inhomogeneities is affected by scat- 
tering. In scattering materials, such as powders, disordered metamaterials or bio- 
logical tissue, multiple scattering on sub-wavelength particles makes light diffuse. 
Recently, we showed that it is possible to construct a wavefront that focuses through 
a solid, strongly scattering object. The focusing wavefront uniquely matches a cer- 
tain configuration of the particles in the medium. To focus light through a turbid 
liquid or living tissue, it is necessary to dynamically adjust the wavefront as the 
particles in the medium move. Here we present three algorithms for constructing a 
wavefront that focuses through a scattering medium. We analyze the dynamic be- 
havior of these algorithms and compare their sensitivity to measurement noise. The 
algorithms are compared both experimentally and using numerical simulations. The 
results are in good agreement with an intuitive model, which may be used to develop 
dynamic diffusion compensators with applications in, for example, light delivery in 
human tissue. 
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1 Introduction 

Materials such as paper, white paint or human tissue are non-transparent be- 
cause of multiple scattering of light [Tf2"f3] . Light propagating in such materials 
is diffuse. Recently, we have shown that coherent light can be focused through 
diffusive media yielding a sharp, intense focus [3]. Starting with the situa- 
tion where a scattering object (a layer of Ti02 pigment with a thickness of 
approximately 20 transport mean free paths) completely destroys the spatial 
coherence of the incident light (Fig. UK, Hfc), we controlled the incident wave- 
front to exactly match scattering in the sample. Afterwards, the transmitted 
light converged to a tight, high contrast focus (Fig. [lb, [Hi). These matched 
wavefronts experience inverse diffusion, that is, they gain spatial coherence by 
travelling through a disordered medium. 

For a given sample of scattering material, there is a unique incident wavefront 
that makes the object optimally focus light to a given point. Like a speckle 
pattern, this wavefront is disordered on the scale of the wavelength of light. 
This wavefront cannot be constructed from a small number of smooth base 
functions, which unfortunately renders the efficient algorithms used in adap- 
tive optics (see e.g. [5]) ineffective. In Ref. pE], we presented an algorithm 
that finds the optimal wavefront when the sample is perfectly stationary and 
the noise level is negligible. To find applications in, for example, fluorescence 
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Fig. 1. Principle and experimental results of inverse wave diffusion, a) A multi- 
ply scattering object destroys the spatial coherence of incident light, b) When the 
same object is illuminated with a specially constructed matching wavefront, the 
transmitted light focuses to a tight spot, c) Recorded intensity transmission of an 
unshaped wave through a 10 /im thick layer of T1O2 pigment, d) Intensity trans- 
mission through the same sample with a shaped wavefront. 



excitation or photodynamic therapy, the wavefront has to be adjusted dy- 
namically as the scatterers in the sample move. In this paper, we present two 
additional algorithms, that dynamically adjust the wavefront to follow changes 
in the sample. The performance of the algorithms is in good agreement with 
numerical simulations and with an analytical model. We show that the new 
algorithms are superior to the original algorithm when the scatterers in the 
sample move or when the initial signal to noise ratio is poor. 

Wave diffusion is a widely encountered physical phenomenon. The use of multi- 
ple scattered waves is the subject of intensive study in the fields of, for instance, 
ultrasound imaging [6fTf8] . radio and microwave antennas [9lfT0] . seismogra- 
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phy [TT] , submarine communication [12] , and surface plasmons [13] . While the 
algorithms discussed in this paper were developed for spatial phase shaping of 
light, they can be used for any type of wave and apply to spatial phase shaping 
as well as to frequency domain phase shaping (also known as coherent control, 
see e.g. [14]) as the concepts are the same. 

This article is organized as follows. First the key concepts of inverse diffusion 
are introduced and the three different algorithms are presented. Then the 
experimental apparatus is explained and the measured typical performance of 
the algorithms is compared. In the subsequent section, we will compare the 
experimental results with numerical simulations and analyze the data in terms 
of noise and stability of the scatterers. Finally, we will analytically explain the 
characteristic features of the different algorithms and discuss their sensitivity 
to noise. 

2 Algorithms for inverse diffusion 

The key elements of an inverse diffusion setup are a multiply scattering sam- 
ple, a spatial light modulator, an optimization algorithm and a detector, as 
shown in Fig. [2J The sample can be anything that scatters light without ab- 
sorbing it. We will consider only samples that are thicker than approximately 
6 transport mean free paths for light. Light transmitted through these samples 
is completely diffuse and the transmitted wavefront is completely scrambled, 
i.e., it has no correlation with the incident wavefront [15J. 

The incident wavefront is constructed using a spatial phase modulator. The 
modulator consists of a 2D-array of pixels that are grouped into N equally 
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sized square segments. A computer sets the phase retardation for each of the 
segments individually to a value between and 2tt. The optimization algorithm 
programs the phase modulator based on the detector output. Since the sample 
completely scrambles the incident wavefront, all segments of the wavefront are 
scattered independently and the optimal wavefront will not be smooth. 

Behind the sample is a detector that provides feedback for the algorithm. The 
detector defines the target area where the intensity is maximized. The field at 
the detector is the result of interference from scattered light originating from 
the different segments of incident wavefront. When the phase of one or more 
segments is changed, the target intensity responds sinusoidally. We sample 
the sine wave by taking 10 measurements. The process of capturing a single 
sine wave and possibly adjusting the phase modulator accordingly is called an 
iteration. 

The amount of control we have over the propagation of light in the disordered 
system is quantified by the signal enhancement. The enhancement rj is defined 

as 



where 1^ is the intensity in the target after optimization and (Iq) is the ensem- 
ble averaged transmitted intensity before optimization. In a perfectly stable 
system, the enhancement is proportional to N [2], meaning that the more 
individual segments are used to shape the incident wavefront, the more light 
is directed to the target. In practice, however, the enhancement is limited by 
the number of iterations that can be performed before the sample changes too 
much. We define the persistence time T p as the decay time of the field autocor- 
relate of the transmitted speckle, which is a measure for the temporal stability 
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Fig. 2. Feedback loop for achieving inverse diffusion. An incident monochromatic 
beam is shaped using a spatial light modulator (1) and projected on a non-transpar- 
ent multiply scattering object (2). A detector (3) detects the amount of transmitted 
light that reaches the target area. A feedback algorithm (4) uses the signal from 
the detector to program the phase modulator. Before the algorithm is started, the 
transmitted light forms a random speckle pattern. The algorithm changes the in- 
cident wave to increase the intensity in the target area. After a few iterations, the 
transmitted light focuses on the target. 

of the sample. The persistence time depends on the type of sample and on en- 
vironmental conditions. Typical values of T p range from a few milliseconds in 
living tissue [16J to hours for solid samples in laboratory conditions. The other 
relevant timescale is the time required for performing a single iteration, Tj. In 
our experiments, we operate the phase modulator at just below 10 Hz and 
take ten measurements for each iteration; we have Tj « 1.2s. 



We will now present three algorithms we used to invert wave diffusion. The 
advantages and disadvantages of the algorithms are discussed briefly and will 
be analyzed in detail later. 
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Fig. 3. Principle used in the three different optimization algorithms, a) For the step- 
wise sequential algorithm, all segments are addressed sequentially (marked squares). 
After the optimal phase is measured for all segments, the modulator is updated to 
construct the optimal wavefront (light gray squares), b) The continuous sequential 
algorithm is equal to the first algorithm, except that the modulator is updated after 
each iteration, c) The partitioning algorithm randomly selects half of the segments 
and adjusts their overall phase. The modulator is updated after each measurement. 



2.0.1 The stepwise sequential algorithm 



The stepwise sequential algorithm that was used in Ref. [I] is very straightfor- 
ward. The computer consecutively cycles the phase of each of the N segments 
from to 2tt. The feedback signal is monitored and the phase for which the 
target intensity is maximal is stored. After all iterations are performed, the 
phase of each segment is set to this optimal value (see Fig. [3^). In absence of 
measurement noise or temporal instability, algorithm 1 is guaranteed to find 
the global maximum in the lowest number of iterations possible. However, 
when NTi ^> T p , the speckle pattern decorrelates before all measurements 
are performed and the algorithm will not work. Therefore, it is important to 
adjust the number of segments to the persistence time. 
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2.0.2 The continuous sequential algorithm 

The continuous sequential algorithm is very similar to the stepwise sequential 
algorithm except for the fact that the phase of each segment is set to its max- 
imum value directly after each measurement (see Fig. [3b). This approach has 
two advantages. First of all, the algorithm runs continuously and dynamically 
follows changes in the sample's scattering behavior. Furthermore, the target 
signal starts to increase directly, which increases the signal to noise ratio of 
successive measurements. It is still necessary to adjust N to the persistence 
time T p . 

2.0.3 The partitioning algorithm 

As an alternative to the two sequential algorithms, we propose a partitioning 
algorithm that requires no a-priori information about the sample's stability. 
Each iteration the phase modulator is divided randomly into two subsets, both 
containing half of the segments (Fig. [3b). The target intensity is maximized 
by changing the phase of one subset with respect to the other. Since the phase 
of half of the segments is changed, the initial increase in intensity will be fast 
and the feedback signal will be maximal. Therefore, this algorithm is expected 
to be less sensitive to noise and to recover from disturbances more rapidly. 

3 Experiment 

The different algorithms were tested experimentally using the setup shown 
in Fig. HI In our case the scattering medium is a 10 /im thick layer of rutile 
TiC>2 pigment jTTj with a mean free path of 0.55 ± 0.1 /im, determined by 
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measuring the total transmission at a wavelength of 632.8 nm. This sample 
is illuminated by a 632.8 nm HeNe laser. The laser beam is expanded and 
spatially modulated by a Holoeye R-2500 liquid crystal light modulator (LCD) 
operating in phase- mostly modulation mode [18]. The shaped beam is focused 
on the sample using a 63x objective with a numerical aperture (NA) of 0.85. A 
20x objective (NA=0.5) images a point that is approximately 3.5 mm behind 
the sample onto a 12-bit CCD camera (Allied Vision Technologies Dolphin 
F-145B). This point is the target area where we want the light to focus. A 
computer integrates the intensity in a circular area with a radius of 20 pixels 
(corresponding to 129 /j,m in the focal plane of the objective). This target 
area is smaller than a typical speckle spot. Using this signal as feedback, the 
computer programs the phase modulator using one of the algorithms described 
above. 

We first run the three different algorithms with iV = 52. Since T p /Ti ^> 52, 
we do not expect to see decoherence effects. In total, 208 iterations were 
performed, which means that the sequential algorithms ran four times con- 
secutively. The results of the optimization procedures is shown in Fig. [5h- 
Although the three algorithms reach the same final enhancement of intensity, 
there are significant differences between the algorithms. The enhancement 
for the stepwise sequential algorithm increases in discrete steps because the 
phase modulator is only reprogrammed every N iterations. During the first 
iV iterations, the target signal is low and the algorithm suffers from noise. 
The saturation enhancement is reached after the second update (after 2N it- 
erations). The continuous sequential algorithm and the partitioning algorithm 
both start updating the wavefront immediately and, therefore, have a higher 
initial increase of the signal. The continuous sequential algorithm is the first 
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Fig. 4. Experimental apparatus used for inverting diffusion. Light from a HeNe laser 
is spatially modulated by a liquid crystal spatial light modulator (SLM). Wave plates 
and a polarizer are used to generate and select the polarization state for which 
the modulator works in phase mostly mode. The shaped beam is focused on the 
sample. A reference detector monitors the total intensity falling on the sample. A 
microscope objective, a polarizer and a CCD-camera are used to detect the intensity 
in the target focus, a few millimeters behind the sample. 



algorithm to reach the saturation enhancement (after A" iterations). The par- 
titioning algorithm has the fastest initial increase in the target signal. It is, 
however, the last algorithm to reach the saturation enhancement since the 
final convergence is very slow. 

When the number of segments in the wavefront is increased, we expect to find a 
higher target intensity. Figure [6^ shows the experimental results for A" = 1804 
on a logarithmic scale. The final intensity enhancement is approximately 40 
times higher than in Fig. [SK A further difference is that in this situation the 
effects of decoherence are no longer negligible. This effect is most clearly visible 
with the stepwise sequential algorithm. The phase modulator is updated after 
each A" iterations and between the updates the intensity decays exponentially 
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Measured enhancement (N=52) 
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Fig. 5. (a) Typical runs of the stepwise sequential algorithm (dotted line), the contin- 
uous sequential algorithm (dashed line) and the partitioning algorithm (solid line). 
All algorithms were run with N = 52. The sequential algorithms were repeated four 
times, (b) Simulation results for N = 52 averaged over 64 runs. The simulation 
captures the main features of the three algorithms, but predicts a higher maximum 
enhancement. 



with a 1/e decay of about T p /Ti = 5000 iterations. 



The convergence behavior of the three algorithms is similar to the experiment 
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shown in Fig. [51 The partitioning algorithm clearly causes a higher signal 
enhancement during the first 1000 iterations. The initial increase in the en- 
hancement is linear with a slope of 0.37. Initially, this linear increase is far 
superior to the quadratic increase obtained with the continuous sequential 
algorithm. 

We conclude that both new algorithms are valuable improvements over the 
original stepwise sequential algorithm. These algorithms are far less sensitive 
to noise and the target signal is kept at a constant value even in the presence 
of decoherence. The partitioning algorithm has the fastest initial increase and 
therefore will recover from disturbances most rapidly. 

4 Simulations 

In order to obtain a better understanding of the effects of noise and fluctu- 
ations on the performance of the different algorithms we perform numerical 
simulations. The disordered medium is represented by a transmission matrix 
with elements drawn from a circular Gaussian distribution (more details on 
the matrix representation can be found below). Decoherence is modelled by 
adding a small perturbation to the transmission matrix after every measure- 
ment. Finally, measurement noise is included by adding a random value to the 
simulated detector signal. 

Figures Eb and[Hb show the simulated enhancement for a system with T p /Tj = 
5000. Every iteration, ten measurements are performed for phase delays be- 
tween and 2tt. To each of these measurements, Gaussian noise with a stan- 
dard deviation of 0.3/ was added. The magnitude of the noise is comparable 
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Fig. 6. (a) Typical runs of stepwise sequential algorithm (dotted line), the continuous 
sequential algorithm (dashed line) and the partitioning algorithm (solid line). All 
algorithms were run with N = 1804. The enhancement are plotted on a logarithmic 
scale, (b) Simulation results for N = 1804 averaged over 64 runs. 



to experimental observations. The three different algorithms were run with 
N = 52 and N = 1804 to simulate the experiments shown in Fig. [5^ and Fig. 

The simulations are in good qualitative agreement with the experimental data. 
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The result for the stepwise sequential algorithm shows that the effects of noise 
and decoherence are simulated realistically. Furthermore, the initial signal in- 
crease and the long time convergence behavior correspond to the measured 
results. The only significant difference is the 20% to 50% higher enhance- 
ment reached in the simulations. A possible explanation for this difference is 
the residual amplitude cross-modulation in our phase modulator. Due to this 
cross-modulation, the optimal wavefront cannot be generated exactly. Further- 
more, the amplitude modulation decreases the accuracy of the measurement 
of the optimal phase. The partitioning algorithm is less sensitive to this last ef- 
fect since the cross-modulation is averaged over many segments with different 
phases. Since the simulations capture the overall behavior of the algorithms 
very well, we can use them to extrapolate to situations with a lot of noise and 
strong decoherence or, on the other hand, to perfectly stable systems. 

5 Analytical expressions for the enhancement 

In this section, we analyze the performance of the algorithms with analytical 
theory and compare these results to the simulations. We describe scattering in 
the sample with the transmission matrix elements, t mn . This matrix couples 
the fields of the incident light and the transmitted light. 

N 

E m = t rnn A n e (2) 

n 

where the <f> n is the phase of the nth segment of the phase modulator. Assuming 
that the modulator is illuminated homogeneously, all incoming channels carry 
the same intensity. We write A n = l/y/N to normalize the total incident 
intensity. Elements Ei, E2, . . . correspond to single scattering channels of the 
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transmitted light. Since we are interested in focusing light to a single spot, 
we need to consider only a single transmission channel, E m . The intensity 
transmitted into channel m is given by 



1 

N 



N 

L 



' L mn z ' 



(3) 



Regardless of the values of the elements t mn of the transmission matrix, the 
intensity \E m \ 2 has its global maximum when the phase modulator exactly 
compensates the phase retardation in the sample for each segment, i.e. <p n = 
— arg(t mn ). The target intensities before optimization (Jo) and after an ideal 
optimization (/ max ) are given by 



1 

N 



N 



(4) 



and 



v 



N 



(5) 



For a disordered medium the elements of t mn are independent and have a 
Gaussian distribution [19.20 .21. 22] . Rewriting Eq. (jSJ) gives 



/ i N i N 

(-^max) ( "77 ^ ] |^mn| |^mfc| "77 ^ ] |^ 
V V n,fc^n iV n 



mn | / > 



7T 



(6) 
(7) 



where the angled brackets denote ensemble averaging over disorder. Eq. (J7J) 
predicts that the expected maximum enhancement for an ideally stable, noise 
free system linearly depends on the number of segments N. For N 3> 1, we 
have rj « tcN/4. 
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5. 1 Performance in fluctuating environments 



In reality, the sample will not be completely stable. Whether this instability 
is due to a drift of the sample position, movement of the scatterers, changing 
humidity or any other cause, the transmission matrix will fluctuate over time. 
In the simulations, we modelled decoherence by repeatedly adding a small 
perturbation to each of the matrix elements. 

tmn > / r r, (tmn + £)' ip) 

yl + <r 

where £ is drawn from a complex Gaussian distribution with mean and 
standard deviation 5. The prefactor normalizes the transformation so that 
(\t\ 2 ) remains constant. By substituting the continuous limit of Eq. (jSJ) in Eq. 
([6]), we find an analytic expression for the effect of decoherence, 



(In) = (Io) 



_ / n \ z 

" '^e-^/Cai)] +0(1) 



(9) 



AN 

where T n is the time that has past since the phase of segment n was measured. 
This simple model explains the exponential decay of the intensity that was 
observed in the measurements (see Figj6k) and predicts a decay time of T p = 
Tt/5 2 . 

We now calculate the maximum enhancement that can be reached with the 
continuous sequential algorithm in the presence of decoherence. Because the 
phases of the segments are measured sequentially, at any given time the values 
for T n are equally spaced between 1 and N. From Eq. we find a maximum 
intensity enhancement of 

v s m = M ^ - 1 +0(1) ' (10) 
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Fig. 7. Theoretical maximum enhancement as a function of coherence time for dif- 
ferent algorithms. The solid lines represent the maximum enhancement that can be 
obtained using the sequential algorithms. The enhancement depends on the number 
of segments used in the algorithm. The dashed line shows the enhancement for the 
partitioning algorithm where N ^> T p /Ti. 

The maximum enhancement for both sequential algorithms is the same. How- 
ever, since the stepwise algorithm only updates the projected wavefront after 
N iterations, the enhancement decreases exponentially between updates. 

In Fig. [7] the enhancement for different values of N is plotted versus T p /Tj. 
When the persistence time is large (T p /Ti ^> N), decoherence effects do not 
play a role and the enhancement linearly depends on N as was seen in Eq. 
([7J. For T p /Ti < N, however, the enhancement decreases because the speckle 
pattern decorrelates before all iterations are performed and the enhancement 
drops to zero. As a consequence, the sequential algorithms only perform opti- 
mal when N is adjusted to T p . When T p is known a-priori, this optimum for 
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N can be found by maximizing Eq. (flOl) . We find 



iVopt = WT p /T h (11) 

where W ~ 2.51 is the solution of exp(W/2) = 1 + W. The maximal enhance- 
ment achievable with sequential algorithms follows by substituting Eq. ffTTT) 
into Eq. (fTOj) and equals r] opt = 0.640T P /Ti. 

With the partitioning algorithm rj increases by 1/2 each iteration of the algo- 
rithm (see appendix [A]). As long as > T p /T^ the enhancement saturates 
at rj — T p /(2Ti) + 1, when the increase is exactly cancelled by the effect of 
decoherence. The most important difference with the sequential algorithms, 
is that the enhancement reached with the partitioning algorithm does not 
depend on N. In Fig. [7J it can be seen that the partitioning algorithm outper- 
forms the sequential algorithms for almost all combinations of T p and N. The 
sequential algorithm only give a slightly higher enhancement when they are 
fine-tuned for a known persistence time (N = 2.51T P /Tj). In most situations, 
T p is not known a-priory or varies over time and the partitioning algorithm 
will be preferable. 

Our analytical results for all three algorithms are supported by numerical sim- 
ulations (see Fig. [HD- The simulations exactly reproduce the theoretical curves 
shown in Fig. [7J For the simulations we used N = 4096, a number that can 
easily be reached with a LCD phase modulator. Again, the partitioning algo- 
rithm can be seen to have good overall performance, whereas the sequential 
algorithms only work well for certain combinations of N and T p . 

In conclusion, the maximum enhancement that can be reached linearly de- 
pends on the sample's persistence time. For the sequential algorithms rj = 
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Fig. 8. Simulated effect of decoherence on the sequential algorithms (solid line) and 
on the partitioning algorithm (dashed line). The simulations are averaged over 25 
runs. Only when TV" ~ 2.51T p /T,, the sequential algorithms perform slightly better 
than the partitioning algorithm. 



0.64T p /Tj, but only when N is precisely adjusted to T p . The partitioning algo- 
rithm has 77 = 0.5T p /Ti, as long as N is large enough. Using these analytical 
relations, the performance of each of the three algorithms in different experi- 
mental situations can easily be estimated. 



6 Effect of Noise 



Measurement noise affects the measured phases. Noise induced errors in the 
phases lead to a reduction of the enhancement, rj. We will now compare the 
signal-to-noise ratio (SNR) of the three different algorithms. In a single iter- 
ation of an algorithm, the phase of one or more segments is varied, while the 
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Table 1 



Signal and noise characteristics of the three algorithms. The rms phase correction 
is a measure for the required sensitivity. 

phase of the other segments is kept constant. The intensity at the detector 
equals 

/($) = I A + I B + 2^lJ^cos ($-$ ), (12) 

where Ib is the intensity at the target originating from the modulated seg- 
ments, I a is the target intensity caused by light coming from the other seg- 
ments, <£> is the phase that is varied and $o is the unknown optimal value for 
the phase. The last term in Eq. ( fl2l) is the signal that is relevant for measuring 
$o- There first two terms constitute a constant bias. 

Table [1] lists the magnitudes of the signal and the bias for each of the three 
algorithms. If the detection system is photon shot noise limited, the noise is 
proportional to the square root of the bias. When, on the other hand, constant 
noise sources such as readout noise or thermal noise are dominant, the SNR 
is directly proportional to the signal magnitude. Table [1] also shows the root 
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mean square (rms) phase correction that is applied during each iteration of 
the corresponding algorithm. The rms phase correction is a measure for the 
required accuracy of the measurements. 

With the stepwise sequential algorithm, I a ~ Io and on average Ib = Iq/N. 
Since the initial diffuse transmission I is low and N can be very high, the 
SNR is low. The continuous sequential algorithm has a higher SNR since the 
overall intensity at the detector increases while the algorithm progresses and 
I A ~ V^o- Assuming the dominant noise source is constant, the SNR will 
increase as the enhancement becomes higher. Therefore, the algorithm can be 
accelerated by decreasing the integration time of the camera as the algorithm 
advances. In case the detection system is photon shot noise limited, the SNR 
remains constant during the optimization since both the signal and the shot 
noise scale as ^/t/Iq. 

The highest SNR is achieved with the partitioning algorithm. Since we always 
change the phase of half of the segments, I a ~ Ib ~ 77/0/2, resulting in a 
maximal signal. Unlike the sequential algorithms, the SNR does not depend 
on N. Therefore the number of segments can be increased without suffering 
from noise. Like with the continuous sequential algorithm, the integration time 
can be adjusted dynamically to optimize the speed/SNR tradeoff. Although 
the partitioning algorithm has the highest SNR, the magnitude of the phase 
corrections decreases as the algorithm progresses. The required accuracy in 
measuring <E> increases at the same pace as the SNR increases. 

The partitioning algorithm is very sensitive to measurement errors since, when 
the measured $ has an error, half of the segments will be programmed with 
the wrong phase. In the extreme case where the error equals 7r the enhancement 
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completely disappears in a single iteration. A simple and effective solution to 
this problem is to keep the previous configuration of the phase modulator 
in memory. When an optimization step causes the signal to decrease, the 
algorithm can revert to the saved configuration. 

7 Conclusion 

Three different algorithms for inverting wave diffusion were presented. The 
algorithms were compared experimentally, with numerical simulations and us- 
ing analytical theory. We found good agreement between experimental data, 
simulations and theory. Moreover, the simulations and theory can be used to 
predict the performance in different experimental situations. 

The effectiveness of the algorithms was quantified by the enhancement. It was 
seen that the enhancement depends on the number of segments iV and the 
relative persistence time T p /Tj. For the sequential algorithms to have optimal 
performance, it is required to adjust iV to match T p . This means that these 
algorithms need a-priori knowledge of the system. The partitioning algorithm 
does not need this knowledge and always performs close to optimal. Moreover, 
the algorithm causes the enhancement to increase the most rapidly of the three 
investigated methods. All in all, this algorithm is a good candidate for applying 
inverse diffusion in instable scattering media such as living tissue. In the future, 
learning algorithms (see e.g. [231211) might be developed to further improve 
the performance of inverse diffusion, for instance by dynamically balancing 
the trade-off between signal to noise ratio and measurement speed. 

The maximum enhancement linearly depends on the number of measurements 
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that can be performed before the speckle pattern decorrelates (T p /Ti). The 
faster the measurements, the higher the enhancement. In our current system, 
the speed is limited by the response time of the LCD. Fast micro mechanical 
phase modulators have a mechanical response time of about 10 /xs (see e.g. 
|25j). which allows a 10 4 times faster operation than with our current system. 
In perfused tissue, a typical decorrelation timescale is 10 ms [IB] , which means 
that an enhancement of about 50 should be possible with currently available 
technology. 
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A Calculation of the performance of the partitioning algorithm 

In this appendix we calculate the development of the enhancement of the 
partitioning algorithm under ideal conditions. During one iteration of the par- 
titioning algorithm, the phase modulator is randomly split into two groups (A 
and B), each containing half of the segments. The relative phase ($) of group 
B is cycled from to 2tc. During this cycle, the target intensity is given by 

/($)= E mA + E mB e^\ (A.l) 
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where E mA is the contribution of the segments in group A to the target field 



E m A = V "^fCmn, (A.2) 



neA 



with 



<^mn — \ / t \ ^"mn& i (A.3) 
V \ J 0/ 

and similar for E m B- The coefficients £ mn are initially random and distributed 
according to a normalized circular Gaussian distribution, meaning that (£) = 
and ((Re£) 2 ) = ((Im£) 2 ) = 1/2. As the algorithm proceeds, the phases n 
are adjusted and the distribution gradually changes to a Rayleigh distribution 
when a high enhancement is reached. The average value (£) increases from 



to \/7r/4 as all contributions are aligned to be in phase. At any moment during 



the optimization, (£) = ^(77 - 1)/(N - 1) and (|£| 2 ) = 1. 

Figure IA.1I gives a graphical representation of a single iteration. Before the 
iteration, E mA and E mB have a different phase. Without loss of generality, we 
choose the phase of {E m A + E m B) to be 0. The intensity before the iteration 
is given by 

-^before 

(Re E mA + Re E mB ) 1 . (A.4) 

After the iteration, $ is set to the value that caused the highest target inten- 
sity, which means that E m & and E m B are now in phase. The target intensity 
then equals 

-^after 

(\E mA \ + \E mB \Y, (A.5) 
which is higher than or equal to before the iteration. 

We now calculate the average intensity gained in a single iteration. We consider 
the regime where already a few iterations have been done (rj ^> 1). In this 
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/V 



b 



Fig. A.l. Complex plane representation of the partitioning algorithm, a) Before the 
iteration the contributions from A and B are not exactly in phase, b) After the 
iteration, the contributions are aligned and the resulting intensity is higher. 



regime, we can approximate 



E mA \ = V( Re E mA? + (Im E mA f (A.6) 

(Im E m A) 2 , . , 

« Re E mA + \ ' . (A.7) 
Ztie ry m A 



Using this result in Eq. (IA.50 gives 



-^after 

(Re E mA + Re E mB \f + (Im E mA f + (Im E mB ) 2 + 

Re E mA 2 Re _E mjB „ , 2 , 

+ 5—= — (ImE mB ) +^-^ — (ImM + 
rte i^ m s -tie i^ m A 

(Im E mA Y (Im E mB )^ (Im £ m Alm ^ m Ji) 2 , . g x 

4(Re E mA ) 2 A(ReE mB f 2Re E mA Re E mB [ ' j 

where the terms on the last line can be neglected. If N ^> 1, 

Re -E m #/Re _E myl « 1. The intensity gain for the iteration AI = / after — /before 

is now found to be 

AI = 2(Im E mA f + 2(Im E mB f (A.9) 

We are primarily interested in the regime N ^> i] ^> 1), where the algorithm 
picks up the main part of the final enhancement. In this regime, (£) <C 1 
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and the probability distribution of £ is still close to the original Gaussian 
distribution. Therefore, ((Im £) 2 ) « 1/2 and it follows from Eq. (1A.2I) that 



Therefore, we expect the intensity enhancement 77 to increase with 1/2 after 
each iteration of the algorithm. With this information, we also calculate the 
typical phase adjustment that is performed in each iteration. From Fig. I A. 1 1 
it follows that the root mean square phase adjustment equals 



When r] approaches its maximum, all contributions are almost completely in 
phase and ((Im £) 2 ) vanishes. In this regime, the algorithm becomes less and 
less effective, as was seen in simulations and experiments (see Fig. [6]). 
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